home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DMLSS.FOR < prev    next >
Text File  |  1985-11-29  |  6KB  |  203 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DMLSS
  5. C
  6. C        PURPOSE
  7. C           SUBROUTINE DMLSS IS THE SECOND STEP IN THE PROCEDURE FOR
  8. C           CALCULATING THE LEAST SQUARES SOLUTION OF MINIMAL LENGTH
  9. C           OF A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH SYMMETRIC
  10. C           POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX.
  11. C
  12. C        USAGE
  13. C           CALL DMLSS(A,N,IRANK,TRAC,INC,RHS,IER)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           A     - COEFFICIENT MATRIX IN FACTORED FORM AS GENERATED
  17. C                   BY SUBROUTINE MFSS FROM INITIALLY GIVEN SYMMETRIC
  18. C                   COEFFICIENT MATRIX A STORED IN N*(N+1)/2 LOCATIONS
  19. C                   A REMAINS UNCHANGED
  20. C                   A MUST BE OF DOUBLE PRECISION
  21. C           N     - DIMENSION OF COEFFICIENT MATRIX
  22. C           IRANK - RANK OF COEFFICIENT MATRIX, CALCULATED BY MEANS OF
  23. C                   SUBROUTINE DMFSS
  24. C           TRAC  - VECTOR OF DIMENSION N CONTAINING THE
  25. C                   SUBSCRIPTS OF PIVOT ROWS AND COLUMNS, I.E. THE
  26. C                   PRODUCT REPRESENTATION IN TRANSPOSITIONS OF THE
  27. C                   PERMUTATION WHICH WAS APPLIED TO ROWS AND COLUMNS
  28. C                   OF A IN THE FACTORIZATION PROCESS
  29. C                   TRAC IS A RESULTANT ARRAY OF SUBROUTINE MFSS
  30. C                   TRAC MUST BE OF DOUBLE PRECISION
  31. C           INC   - INPUT VARIABLE WHICH SHOULD CONTAIN THE VALUE ZERO
  32. C                   IF THE SYSTEM OF SIMULTANEOUS EQUATIONS IS KNOWN
  33. C                   TO BE COMPATIBLE AND A NONZERO VALUE OTHERWISE
  34. C           RHS   - VECTOR OF DIMENSION N CONTAINING THE RIGHT HAND SIDE
  35. C                   ON RETURN RHS CONTAINS THE MINIMAL LENGTH SOLUTION
  36. C                   RHS MUST BE OF DOUBLE PRECISION
  37. C           IER   - RESULTANT ERROR PARAMETER
  38. C                   IER = 0 MEANS NO ERRORS
  39. C                   IER =-1 MEANS N AND/OR IRANK IS NOT POSITIVE AND/OR
  40. C                           IRANK IS GREATER THAN N
  41. C                   IER = 1 MEANS THE FACTORIZATION CONTAINED IN A HAS
  42. C                           ZERO DIVISORS AND/OR TRAC CONTAINS
  43. C                           VALUES OUTSIDE THE FEASIBLE RANGE 1 UP TO N
  44. C
  45. C        REMARKS
  46. C           THE MINIMAL LENGTH SOLUTION IS PRODUCED IN THE STORAGE
  47. C           LOCATIONS OCCUPIED BY THE RIGHT HAND SIDE.
  48. C           SUBROUTINE DMLSS DOES TAKE CARE OF THE PERMUTATION
  49. C           WHICH WAS APPLIED TO ROWS AND COLUMNS OF A.
  50. C           OPERATION IS BYPASSED IN CASE OF A NON POSITIVE VALUE
  51. C           OF IRANK
  52. C
  53. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  54. C           NONE
  55. C
  56. C        METHOD
  57. C           LET T, U, TU BE THE COMPONENTS OF THE FACTORIZATION OF A,
  58. C           AND LET THE RIGHT HAND SIDE BE PARTITIONED INTO A FIRST
  59. C           PART X1 OF DIMENSION IRANK AND A SECOND PART X2 OF DIMENSION
  60. C           N-IRANK. THEN THE FOLLOWING OPERATIONS ARE APPLIED IN
  61. C           SEQUENCE
  62. C           (1) INTERCHANGE RIGHT HAND SIDE
  63. C           (2) X1 = X1 + U * X2
  64. C           (3) X2 =-TRANSPOSE(U) * X1
  65. C           (4) X2 = INVERSE(TU) * INVERSE(TRANSPOSE(TU)) * X2
  66. C           (5) X1 = X1 + U * X2
  67. C           (6) X1 = INVERSE(T) * INVERSE(TRANSPOSE(T)) * X1
  68. C           (7) X2 =-TRANSPOSE(U) * X1
  69. C           (8) X2 = INVERSE(TU) * INVERSE(TRANSPOSE(TU)) * X2
  70. C           (9) X1 = X1 + U * X2
  71. C           (10)X2 = TRANSPOSE(U) * X1
  72. C           (11) REINTERCHANGE CALCULATED SOLUTION
  73. C           IF THE SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS IS SPECIFIED
  74. C           TO BE COMPATIBLE THEN STEPS (2), (3), (4) AND (5) ARE
  75. C           CANCELLED.
  76. C           IF THE COEFFICIENT MATRIX HAS RANK N, THEN THE ONLY STEPS
  77. C           PERFORMED ARE (1), (6) AND (11).
  78. C
  79. C     ..................................................................
  80. C
  81.       SUBROUTINE DMLSS(A,N,IRANK,TRAC,INC,RHS,IER)
  82. C
  83. C
  84. C        DIMENSIONED DUMMY VARIABLES
  85.       DIMENSION A(1),TRAC(1),RHS(1)
  86.       DOUBLE PRECISION SUM,A,RHS,TRAC,HOLD
  87. C
  88. C        TEST OF SPECIFIED DIMENSIONS
  89.       IDEF=N-IRANK
  90.       IF(N)33,33,1
  91.     1 IF(IRANK)33,33,2
  92.     2 IF(IDEF)33,3,3
  93. C
  94. C        CALCULATE AUXILIARY VALUES
  95.     3 ITE=IRANK*(IRANK+1)/2
  96.       IX2=IRANK+1
  97.       NP1=N+1
  98.       IER=0
  99. C
  100. C        INTERCHANGE RIGHT HAND SIDE
  101.       JJ=1
  102.       II=1
  103.     4 DO 6 I=1,N
  104.       J=TRAC(II)
  105.       IF(J)31,31,5
  106.     5 HOLD=RHS(II)
  107.       RHS(II)=RHS(J)
  108.       RHS(J)=HOLD
  109.     6 II=II+JJ
  110.       IF(JJ)32,7,7
  111. C
  112. C        PERFORM STEP 2 IF NECESSARY
  113.     7 ISW=1
  114.       IF(INC*IDEF)8,28,8
  115. C
  116. C        CALCULATE X1 = X1 + U * X2
  117.     8 ISTA=ITE
  118.       DO 10 I=1,IRANK
  119.       ISTA=ISTA+1
  120.       JJ=ISTA
  121.       SUM=0.D0
  122.       DO 9 J=IX2,N
  123.       SUM=SUM+A(JJ)*RHS(J)
  124.     9 JJ=JJ+J
  125.    10 RHS(I)=RHS(I)+SUM
  126.       GOTO(11,28,11),ISW
  127. C
  128. C        CALCULATE X2 = TRANSPOSE(U) * X1
  129.    11 ISTA=ITE
  130.       DO 15 I=IX2,N
  131.       JJ=ISTA
  132.       SUM=0.D0
  133.       DO 12 J=1,IRANK
  134.       JJ=JJ+1
  135.    12 SUM=SUM+A(JJ)*RHS(J)
  136.       GOTO(13,13,14),ISW
  137.    13 SUM=-SUM
  138.    14 RHS(I)=SUM
  139.    15 ISTA=ISTA+I
  140.       GOTO(16,29,30),ISW
  141. C
  142. C        INITIALIZE STEP (4) OR STEP (8)
  143.    16 ISTA=IX2
  144.       IEND=N
  145.       JJ=ITE+ISTA
  146. C
  147. C        DIVISION OF X1 BY TRANSPOSE OF TRIANGULAR MATRIX
  148.    17 SUM=0.D0
  149.       DO 20 I=ISTA,IEND
  150.       IF(A(JJ))18,31,18
  151.    18 RHS(I)=(RHS(I)-SUM)/A(JJ)
  152.       IF(I-IEND)19,21,21
  153.    19 JJ=JJ+ISTA
  154.       SUM=0.D0
  155.       DO 20 J=ISTA,I
  156.       SUM=SUM+A(JJ)*RHS(J)
  157.    20 JJ=JJ+1
  158. C
  159. C        DIVISION OF X1 BY TRIANGULAR MATRIX
  160.    21 SUM=0.D0
  161.       II=IEND
  162.       DO 24 I=ISTA,IEND
  163.       RHS(II)=(RHS(II)-SUM)/A(JJ)
  164.       IF(II-ISTA)25,25,22
  165.    22 KK=JJ-1
  166.       SUM=0.D0
  167.       DO 23 J=II,IEND
  168.       SUM=SUM+A(KK)*RHS(J)
  169.    23 KK=KK+J
  170.       JJ=JJ-II
  171.    24 II=II-1
  172.    25 IF(IDEF)26,30,26
  173.    26 GOTO(27,11,8),ISW
  174. C
  175. C        PERFORM STEP (5)
  176.    27 ISW=2
  177.       GOTO 8
  178. C
  179. C        PERFORM STEP (6)
  180.    28 ISTA=1
  181.       IEND=IRANK
  182.       JJ=1
  183.       ISW=2
  184.       GOTO 17
  185. C
  186. C        PERFORM STEP (8)
  187.    29 ISW=3
  188.       GOTO 16
  189. C
  190. C        REINTERCHANGE CALCULATED SOLUTION
  191.    30 II=N
  192.       JJ=-1
  193.       GOTO 4
  194. C
  195. C        ERROR RETURN IN CASE OF ZERO DIVISOR
  196.    31 IER=1
  197.    32 RETURN
  198. C
  199. C        ERROR RETURN IN CASE OF ILLEGAL DIMENSION
  200.    33 IER=-1
  201.       RETURN
  202.       END
  203.